{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"id": "0d79d906",
"metadata": {
"nbsphinx": "hidden",
"execution": {
"iopub.execute_input": "2026-09-29T20:07:05.770087Z",
"iopub.status.busy": "2026-09-29T20:07:05.769981Z",
"iopub.status.idle": "2026-09-29T20:07:06.409022Z",
"shell.execute_reply": "2026-09-29T20:07:06.408100Z"
}
},
"outputs": [],
"source": [
"%matplotlib inline\n",
"%config InlineBackend.figure_format = 'svg'\n",
"\n",
"import warnings \n",
"warnings.filterwarnings(\"ignore\") #ignore some matplotlib warnings\n",
"\n",
"import numpy as np\n",
"\n",
"from mpi4py import MPI\n",
"\n",
"from triqs.plot.mpl_interface import plt, oplot, oplotr, oploti\n",
"plt.rcParams[\"figure.figsize\"] = (6,4) # set default size for all figures"
]
},
{
"cell_type": "markdown",
"id": "e10b090f",
"metadata": {},
"source": [
"# One-boson retarded interaction\n",
"\n",
"This tutorial shows how to include a retarded (dynamical) interaction in a `triqs_xca` calculation. It is adapted from the [TRIQS/ctseg tutorial of the same name](https://triqs.github.io/ctseg/latest/tutorials/One-boson%20retarded%20interaction.html), which solves the model below with continuous-time quantum Monte Carlo, so that the two solvers can be compared directly.\n",
"\n",
"The model is the smallest one with a retarded interaction: a single fermionic level (0) hybridized with one bath level (1) and coupled linearly to one bosonic mode,\n",
"\n",
"$$\n",
"H = \\epsilon_0 c^\\dagger_0 c_0 + \\epsilon_1 c^\\dagger_1 c_1 \n",
"+ V (c^\\dagger_1 c_0 + c^\\dagger_0 c_1)\n",
"+ g (b + b^\\dagger) c^\\dagger_0 c_0\n",
"+ \\omega_0 b^\\dagger b\n",
"\\, .\n",
"$$\n",
"\n",
"Integrating out the bath level and the boson leaves an action for the level 0 alone,\n",
"\n",
"$$\n",
"\\mathcal{S} = \n",
"\\int_0^\\beta d\\tau \\, H_{0}(\\tau) \n",
"+ \n",
"\\iint_0^\\beta d\\tau d\\tau' \\, c_0^\\dagger(\\tau) \\Delta(\\tau - \\tau') c_0(\\tau')\n",
"+\n",
"\\frac{1}{2} \\iint_0^\\beta d\\tau d\\tau' \\, n_0(\\tau) \\, \\mathcal{U}(\\tau - \\tau') \\, n_0(\\tau')\n",
"\\, ,\n",
"$$\n",
"\n",
"with $H_{0} = \\epsilon_0 c^\\dagger_0 c_0$, $n_0 = c^\\dagger_0 c_0$, the hybridization function \n",
"$\\Delta(i\\nu_n) = \\frac{V^2}{i\\nu_n - \\epsilon_1}$,\n",
"and the retarded interaction $\\mathcal{U}(i\\omega_n) = - 2 g^2 \\frac{\\omega_0}{\\omega_0^2 - (i\\omega_n)^2} \\,.$\n",
"\n",
"`triqs_xca` expands in the hybridization and in the retarded interaction on the same footing: the expansion order counts hybridization lines and interaction lines together in the pseudo-particle self-energy diagrams. Below we compare the results of the first four expansion orders (the first two being the non-crossing and one-crossing approximations, NCA and OCA) for the single particle Green's function $G_0(\\tau) = - \\langle \\mathcal{T} c_0(\\tau) c_0^\\dagger \\rangle$ and the density correlator $\\chi(\\tau) = \\langle \\mathcal{T} n_0(\\tau) n_0 \\rangle$ with exact diagonalization of $H$. The bosonic Hilbert space is infinite, so the diagonalization works in a truncated boson number space and is converged with respect to the truncation."
]
},
{
"cell_type": "markdown",
"id": "fece0a6c",
"metadata": {},
"source": [
"## Solving the retarded action with `triqs_xca`\n",
"\n",
"**1.** We start from the parameters of the model. The bosonic mode frequency $\\omega_0 = 1$ is the unit of energy and the coupling $g = 1/2$ is in the intermediate coupling regime (at the stronger coupling $g = 1$ of the ctseg tutorial the order-by-order convergence is slow). The chemical potential $\\mu = -g^2/\\omega_0$ compensates the polaron shift of the level, so that $\\epsilon_0 = 0$ corresponds to half filling."
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "20d251ed",
"metadata": {
"execution": {
"iopub.execute_input": "2026-09-29T20:07:06.412182Z",
"iopub.status.busy": "2026-09-29T20:07:06.411996Z",
"iopub.status.idle": "2026-09-29T20:07:06.415330Z",
"shell.execute_reply": "2026-09-29T20:07:06.413962Z"
}
},
"outputs": [],
"source": [
"eps0, eps1, V, g, omega0 = -0.1, +0.1, 0.25, 0.5, 1.\n",
"mu = -g**2 / omega0 # For half-filling at eps0 = 0\n",
"beta = 10."
]
},
{
"cell_type": "markdown",
"id": "80378dfb",
"metadata": {},
"source": [
"\n",
"\n",
"**2.** Next we construct a solver for the local Hamiltonian $H_0 - \\mu n_0$. The imaginary time response functions are represented in the Discrete Lehmann Representation (DLR), with accuracy `eps` and frequency cut-off `w_max` chosen such that the spectra of $\\Delta$, $\\mathcal{U}$ and $G_0$ are covered."
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "9331251c",
"metadata": {
"execution": {
"iopub.execute_input": "2026-09-29T20:07:06.418266Z",
"iopub.status.busy": "2026-09-29T20:07:06.418160Z",
"iopub.status.idle": "2026-09-29T20:07:07.642653Z",
"shell.execute_reply": "2026-09-29T20:07:07.641527Z"
}
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"Starting run with 8 MPI rank(s) at : 2026-09-29 22:07:06.496281\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"____ ____________ _____\n",
"\\ \\/ /\\_ ___ \\ / _ \\\n",
" \\ / / \\ \\/ / /_\\ \\\n",
" / \\ \\ \\____/ | \\\n",
"/___/\\ \\ \\______ /\\____|__ /\n",
" \\_/ \\/ \\/ [github.com/TRIQS/xca]\n",
"\n",
"beta = 10.0, w_max = 5.0, eps = 1e-10, N_DLR = 24\n",
"AtomDiagReal: dim 2 with 2 subspaces dims [1] freq [2] E min/max +0.00E+00/+1.50E-01\n"
]
}
],
"source": [
"from triqs.operators import n\n",
"from triqs_xca import BlockSparseSolver as Solver\n",
"\n",
"S = Solver(\n",
" H_loc=(eps0 - mu) * n('0', 0),\n",
" beta=beta, # Inverse temperature\n",
" gf_struct=[('0', 1)], # Green's function structure 1st index: name, 2nd index: dimension of subspace\n",
" eps=1e-10, # Accuracy of Discrete Lehmann Representation (DLR) used for imaginary time response functions\n",
" w_max=5.0, # DLR frequency cut-off (the spectrum of the model must be in the range [-w_max, +w_max])\n",
" )"
]
},
{
"cell_type": "markdown",
"id": "fdf6b3bc",
"metadata": {},
"source": [
"\n",
"\n",
"**3.** The hybridization function is stored in `S.Delta_tau` on the DLR imaginary time mesh of the solver. We set it from its Matsubara frequency expression $\\Delta(i\\nu_n) = \\frac{V^2}{i\\nu_n - \\epsilon_1}$,"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "ad60cd2e",
"metadata": {
"execution": {
"iopub.execute_input": "2026-09-29T20:07:07.645633Z",
"iopub.status.busy": "2026-09-29T20:07:07.645271Z",
"iopub.status.idle": "2026-09-29T20:07:07.649775Z",
"shell.execute_reply": "2026-09-29T20:07:07.648354Z"
}
},
"outputs": [],
"source": [
"from triqs.gfs import make_gf_dlr_imfreq, make_gf_dlr_imtime, inverse, iOmega_n\n",
"\n",
"Delta_iw = make_gf_dlr_imfreq(S.Delta_tau['0'])\n",
"Delta_iw << V**2 * inverse(iOmega_n - eps1)\n",
"S.Delta_tau['0'] << make_gf_dlr_imtime(Delta_iw);"
]
},
{
"cell_type": "markdown",
"id": "1d5073d8",
"metadata": {},
"source": [
"\n",
"\n",
"**4.** The retarded interaction is a bosonic function of imaginary time, so we construct it on a bosonic DLR mesh with the same $\\beta$, `eps` and `w_max` as the solver, again from its Matsubara frequency expression $\\mathcal{U}(i\\omega_n) = - 2 g^2 \\frac{\\omega_0}{\\omega_0^2 - (i\\omega_n)^2}$.\n",
"\n",
"Retarded interactions are passed to the solver with `set_dynamic_interactions`, as a list of operators $O_i$ together with a matrix valued kernel $D_{ij}(\\tau)$, which the solver interprets as the action term $\\iint_0^\\beta d\\tau d\\tau' \\sum_{ij} O_i^\\dagger(\\tau) \\, D_{ij}(\\tau - \\tau') \\, O_j(\\tau')$ without a factor $\\frac{1}{2}$. Comparing with the action $\\mathcal{S}$ above, the single operator $O = n_0$ therefore comes with the kernel $D(\\tau) = \\frac{1}{2} \\mathcal{U}(\\tau)$."
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "a1e2a497",
"metadata": {
"execution": {
"iopub.execute_input": "2026-09-29T20:07:07.652061Z",
"iopub.status.busy": "2026-09-29T20:07:07.651937Z",
"iopub.status.idle": "2026-09-29T20:07:07.661270Z",
"shell.execute_reply": "2026-09-29T20:07:07.660111Z"
}
},
"outputs": [],
"source": [
"from triqs.gfs import Gf, MeshDLRImTime\n",
"\n",
"b_mesh = MeshDLRImTime(beta=beta, statistic='Boson', eps=S.mesh_tau.eps, w_max=S.mesh_tau.w_max, symmetrize=False)\n",
"\n",
"U_tau = Gf(mesh=b_mesh, target_shape=[1, 1])\n",
"U_iw = make_gf_dlr_imfreq(U_tau)\n",
"U_iw << -2 * g**2 * omega0 * inverse(omega0**2 - iOmega_n*iOmega_n)\n",
"U_tau << make_gf_dlr_imtime(U_iw)\n",
"\n",
"S.set_dynamic_interactions(dynint_ops=[n('0', 0)], dynint_tau=0.5 * U_tau)"
]
},
{
"cell_type": "markdown",
"id": "f76d1e79",
"metadata": {},
"source": [
"\n",
"\n",
"**5.** Finally we run the solver for expansion orders one to four, and evaluate the density correlator $\\chi(\\tau) = \\langle \\mathcal{T} n_0(\\tau) n_0 \\rangle$ at the same order with `eval_one_time_correlator`. For plotting, the DLR results are interpolated onto a dense imaginary time mesh.\n",
"\n",
"The number of diagrams grows rapidly with the order, and the fourth order dominates the run time by far. The diagram evaluation is MPI parallel, and this notebook was executed with eight MPI ranks."
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "decdf2e8",
"metadata": {
"execution": {
"iopub.execute_input": "2026-09-29T20:07:07.663924Z",
"iopub.status.busy": "2026-09-29T20:07:07.663792Z",
"iopub.status.idle": "2026-09-29T20:07:12.630371Z",
"shell.execute_reply": "2026-09-29T20:07:12.629294Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Hybridization: using DLR expansion with N_poles = 24\n",
"iter = 1, diff_G = 5.16E-01, Z-1 = -1.88E-13, eta = 1.59E-01\n",
"iter = 2, diff_G = 9.77E-02, Z-1 = -3.69E-13, eta = 1.97E-01\n",
"iter = 3, diff_G = 1.43E-02, Z-1 = -5.22E-14, eta = 2.02E-01\n",
"iter = 4, diff_G = 1.56E-03, Z-1 = -5.33E-15, eta = 2.03E-01\n",
"iter = 5, diff_G = 9.97E-05, Z-1 = -5.55E-16, eta = 2.03E-01\n",
"iter = 6, diff_G = 4.28E-06, Z-1 = +2.22E-16, eta = 2.03E-01\n",
"iter = 7, diff_G = 1.33E-07, Z-1 = +0.00E+00, eta = 2.03E-01\n",
"Converged after 7 iterations with diff_G = 1.33E-07 < tol = 1.00E-06\n",
"\n",
"Timing: incl. excl.\n",
"----------------------------------------------------\n",
"Adapol hybridization fit: 0.000 0.000 0.3% |\n",
"AtomDiag Init: 0.001 0.001 2.0% ||\n",
"DiagramEvaluator Init: 0.001 0.001 1.6% ||\n",
"Dyson: 0.008 0.001 1.7% ||\n",
" Setup: 0.005 0.005 10.5% |---|\n",
" Solve: 0.002 0.002 4.7% |-|\n",
"Sigma: 0.005 0.000 0.4% |\n",
" Order 1: 0.005 0.005 10.7% |---|\n",
"Single-particle Gf: 0.000 0.000 0.0% |\n",
" Order 1: 0.000 0.000 0.3% |\n",
"Other: 0.032 0.032 67.6% |--------------------------|\n",
"----------------------------------------------------\n",
"Total: 0.048 100.0%\n",
"\n",
"Adapol: compressing 24 poles, tol 1.00E-07 on the imaginary time L2 error\n",
" beta = 10, 24 Matsubara points, linear residue fit\n",
"Adapol: [1/2 bracket] growing the AAA step count until the tolerance is met\n",
" 2 steps -> 3 poles, error 2.28E-11 <= tol ok AAA residual 2.48E-11\n",
"Adapol: [2/2 bisect] smallest step count still meeting the tolerance, candidates 1..2\n",
" 1 steps -> 1 poles, error 3.06E-02 > tol fail AAA residual 1.01E-01 candidates 2..2\n",
"Adapol: done -- 3 poles (2 AAA steps), error 2.28E-11 <= tol 1.00E-07, 2 pipeline passes, 2 AAA runs\n",
"Adapol: Fit error = 2.28E-11 < tol = 1.00E-07, N_poles = 3\n",
"iter = 1, diff_G = 2.70E-02, Z-1 = -1.15E-13, eta = 2.11E-01\n",
"iter = 2, diff_G = 8.20E-03, Z-1 = -3.75E-14, eta = 2.13E-01\n",
"iter = 3, diff_G = 6.05E-04, Z-1 = -3.33E-15, eta = 2.13E-01\n",
"iter = 4, diff_G = 5.84E-05, Z-1 = -2.22E-16, eta = 2.13E-01\n",
"iter = 5, diff_G = 2.68E-06, Z-1 = +4.44E-16, eta = 2.13E-01\n",
"iter = 6, diff_G = 9.30E-08, Z-1 = +4.44E-16, eta = 2.13E-01\n",
"Converged after 6 iterations with diff_G = 9.30E-08 < tol = 1.00E-06\n",
"\n",
"Timing: incl. excl.\n",
"----------------------------------------------------\n",
"Adapol hybridization fit: 0.017 0.017 23.4% |--------|\n",
"AtomDiag Init: 0.001 0.001 1.3% ||\n",
"DiagramEvaluator Init: 0.001 0.001 1.1% |\n",
"Dyson: 0.009 0.001 1.6% ||\n",
" Setup: 0.005 0.005 7.0% |--|\n",
" Solve: 0.003 0.003 4.5% |-|\n",
"One-time correlator: 0.000 0.000 0.0% |\n",
" Order 1: 0.000 0.000 0.1% |\n",
"Sigma: 0.010 0.000 0.4% |\n",
" Order 1: 0.006 0.006 8.9% |---|\n",
" Order 2: 0.003 0.003 4.0% |-|\n",
"Single-particle Gf: 0.001 0.000 0.0% |\n",
" Order 1: 0.000 0.000 0.3% |\n",
" Order 2: 0.000 0.000 0.5% |\n",
"Other: 0.034 0.034 46.7% |------------------|\n",
"----------------------------------------------------\n",
"Total: 0.072 100.0%\n",
"\n",
"Adapol: compressing 24 poles, tol 1.00E-07 on the imaginary time L2 error\n",
" beta = 10, 24 Matsubara points, linear residue fit\n",
"Adapol: [1/2 bracket] growing the AAA step count until the tolerance is met\n",
" 2 steps -> 3 poles, error 2.28E-11 <= tol ok AAA residual 2.48E-11\n",
"Adapol: [2/2 bisect] smallest step count still meeting the tolerance, candidates 1..2\n",
" 1 steps -> 1 poles, error 3.06E-02 > tol fail AAA residual 1.01E-01 candidates 2..2\n",
"Adapol: done -- 3 poles (2 AAA steps), error 2.28E-11 <= tol 1.00E-07, 2 pipeline passes, 2 AAA runs\n",
"Adapol: Fit error = 2.28E-11 < tol = 1.00E-07, N_poles = 3\n",
"iter = 1, diff_G = 6.35E-03, Z-1 = -2.58E-14, eta = 2.14E-01\n",
"iter = 2, diff_G = 1.39E-03, Z-1 = -7.77E-15, eta = 2.15E-01\n",
"iter = 3, diff_G = 7.70E-05, Z-1 = +0.00E+00, eta = 2.15E-01\n",
"iter = 4, diff_G = 6.39E-06, Z-1 = -4.44E-16, eta = 2.15E-01\n",
"iter = 5, diff_G = 2.37E-07, Z-1 = +0.00E+00, eta = 2.15E-01\n",
"Converged after 5 iterations with diff_G = 2.37E-07 < tol = 1.00E-06\n",
"\n",
"Timing: incl. excl.\n",
"----------------------------------------------------\n",
"Adapol hybridization fit: 0.019 0.019 10.8% |---|\n",
"AtomDiag Init: 0.001 0.001 0.5% |\n",
"DiagramEvaluator Init: 0.001 0.001 0.5% |\n",
"Dyson: 0.011 0.001 0.8% |\n",
" Setup: 0.005 0.005 2.8% ||\n",
" Solve: 0.004 0.004 2.3% ||\n",
"One-time correlator: 0.001 0.000 0.0% |\n",
" Order 1: 0.000 0.000 0.1% |\n",
" Order 2: 0.000 0.000 0.2% |\n",
"Sigma: 0.098 0.000 0.2% |\n",
" Order 1: 0.008 0.008 4.5% |-|\n",
" Order 2: 0.006 0.006 3.3% ||\n",
" Order 3: 0.084 0.084 47.0% |------------------|\n",
"Single-particle Gf: 0.014 0.000 0.0% |\n",
" Order 1: 0.000 0.000 0.2% |\n",
" Order 2: 0.001 0.001 0.4% |\n",
" Order 3: 0.012 0.012 7.0% |--|\n",
"Other: 0.035 0.035 19.3% |-------|\n",
"----------------------------------------------------\n",
"Total: 0.179 100.0%\n",
"\n",
"Adapol: compressing 24 poles, tol 1.00E-07 on the imaginary time L2 error\n",
" beta = 10, 24 Matsubara points, linear residue fit\n",
"Adapol: [1/2 bracket] growing the AAA step count until the tolerance is met\n",
" 2 steps -> 3 poles, error 2.28E-11 <= tol ok AAA residual 2.48E-11\n",
"Adapol: [2/2 bisect] smallest step count still meeting the tolerance, candidates 1..2\n",
" 1 steps -> 1 poles, error 3.06E-02 > tol fail AAA residual 1.01E-01 candidates 2..2\n",
"Adapol: done -- 3 poles (2 AAA steps), error 2.28E-11 <= tol 1.00E-07, 2 pipeline passes, 2 AAA runs\n",
"Adapol: Fit error = 2.28E-11 < tol = 1.00E-07, N_poles = 3\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"iter = 1, diff_G = 1.15E-03, Z-1 = -8.99E-15, eta = 2.15E-01\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"iter = 2, diff_G = 3.05E-04, Z-1 = -1.55E-15, eta = 2.15E-01\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"iter = 3, diff_G = 1.39E-05, Z-1 = -4.44E-16, eta = 2.15E-01\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"iter = 4, diff_G = 8.82E-07, Z-1 = -3.33E-16, eta = 2.15E-01\n",
"Converged after 4 iterations with diff_G = 8.82E-07 < tol = 1.00E-06\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
"Timing: incl. excl.\n",
"----------------------------------------------------\n",
"Adapol hybridization fit: 0.022 0.022 0.5% |\n",
"AtomDiag Init: 0.001 0.001 0.0% |\n",
"DiagramEvaluator Init: 0.001 0.001 0.0% |\n",
"Dyson: 0.012 0.002 0.0% |\n",
" Setup: 0.005 0.005 0.1% |\n",
" Solve: 0.005 0.005 0.1% |\n",
"One-time correlator: 0.013 0.000 0.0% |\n",
" Order 1: 0.000 0.000 0.0% |\n",
" Order 2: 0.001 0.001 0.0% |\n",
" Order 3: 0.012 0.012 0.3% |\n",
"Sigma: 3.009 0.001 0.0% |\n",
" Order 1: 0.009 0.009 0.2% |\n",
" Order 2: 0.008 0.008 0.2% |\n",
" Order 3: 0.150 0.150 3.7% ||\n",
" Order 4: 2.841 2.841 70.0% |---------------------------|\n",
"Single-particle Gf: 0.965 0.000 0.0% |\n",
" Order 1: 0.001 0.001 0.0% |\n",
" Order 2: 0.001 0.001 0.0% |\n",
" Order 3: 0.026 0.026 0.6% |\n",
" Order 4: 0.938 0.938 23.1% |--------|\n",
"Other: 0.036 0.036 0.9% |\n",
"----------------------------------------------------\n",
"Total: 4.059 100.0%\n",
"\n"
]
}
],
"source": [
"from triqs.gfs import make_gf_imtime\n",
"\n",
"n_tau = 401\n",
"results = {}\n",
"\n",
"for order in [1, 2, 3, 4]:\n",
" S.solve(max_order=order, tol=1e-6, maxiter=50)\n",
" chi_tau = -S.eval_one_time_correlator(S.G, max_order=order, ops_tau=[n('0', 0)], ops_0=[n('0', 0)]) # returns -\n",
" results[order] = dict(\n",
" G_tau=make_gf_imtime(S.G_tau['0'], n_tau=n_tau), \n",
" chi_tau=make_gf_imtime(chi_tau, n_tau=n_tau))"
]
},
{
"cell_type": "markdown",
"id": "5f7d2400",
"metadata": {},
"source": [
"## Reference solution from exact diagonalization (using **pyED**)\n",
"\n",
"For the reference solution we diagonalize the full Hamiltonian $H$ with two fermions and one boson (truncated to `Nb_max` boson quanta, which is fully converged for the parameters used here) using [pyED](https://github.com/HugoStrand/pyed), and evaluate $G_0(\\tau)$ and $\\chi(\\tau)$ on the same dense imaginary time mesh."
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "6a561c2a",
"metadata": {
"execution": {
"iopub.execute_input": "2026-09-29T20:07:12.633354Z",
"iopub.status.busy": "2026-09-29T20:07:12.633211Z",
"iopub.status.idle": "2026-09-29T20:07:12.651029Z",
"shell.execute_reply": "2026-09-29T20:07:12.650110Z"
}
},
"outputs": [],
"source": [
"from pyed.SparseExactDiagonalization import SparseExactDiagonalization\n",
"from pyed.SparseMatrixFockStates import SparseMatrixFermiBoseCreationOperators\n",
"\n",
"def get_ed_ref(eps0, eps1, V, g, omega0, mesh_tau, Nb_max=10):\n",
"\n",
" ops = SparseMatrixFermiBoseCreationOperators(Nf=2, Nb=1, Nb_max=Nb_max)\n",
" \n",
" c0, c1, b = ops.c_dag[0].getH(), ops.c_dag[1].getH(), ops.b_dag[0].getH()\n",
" n0, n1, nb = c0.getH() * c0, c1.getH() * c1, b.getH() * b\n",
"\n",
" H = eps0 * n0 + eps1 * n1 + V * (c1.getH() * c0 + c0.getH() * c1) + g * (b + b.getH()) * n0 + omega0 * nb\n",
"\n",
" ed = SparseExactDiagonalization(H, mesh_tau.beta)\n",
" tau = np.array([float(t) for t in mesh_tau])\n",
"\n",
" g_tau = Gf(mesh=mesh_tau, target_shape=[1, 1])\n",
" g_tau.data[:, 0, 0] = ed.get_tau_greens_function_component(tau, c0, c0.getH())\n",
"\n",
" chi_tau = Gf(mesh=mesh_tau, target_shape=[1, 1])\n",
" chi_tau.data[:, 0, 0] = -ed.get_tau_greens_function_component(tau, n0, n0) # pyed returns -\n",
" \n",
" return g_tau, chi_tau\n",
"\n",
"g_tau_ed, chi_tau_ed = get_ed_ref(eps0 - mu, eps1, V, g, omega0, results[1]['G_tau'].mesh)"
]
},
{
"cell_type": "markdown",
"id": "96882020",
"metadata": {},
"source": [
"## Single particle Green's function\n",
"\n",
"The first order (NCA) result captures the shape of $G_0(\\tau)$, and every additional order moves the result closer to the exact diagonalization reference."
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "232adc7c",
"metadata": {
"execution": {
"iopub.execute_input": "2026-09-29T20:07:12.654458Z",
"iopub.status.busy": "2026-09-29T20:07:12.654265Z",
"iopub.status.idle": "2026-09-29T20:07:12.826714Z",
"shell.execute_reply": "2026-09-29T20:07:12.825728Z"
}
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"def label(order):\n",
" return f'XCA order {order}' + {1: ' (NCA)', 2: ' (OCA)'}.get(order, '')\n",
"\n",
"for order, r in results.items():\n",
" oplotr(r['G_tau'], label=label(order))\n",
"oplotr(g_tau_ed, 'k--', label='pyED')\n",
"\n",
"plt.xlabel(r'$\\tau$'); plt.ylabel(r'$G_0(\\tau)$'); plt.ylim(top=0); plt.legend(loc='best');"
]
},
{
"cell_type": "markdown",
"id": "d1006365",
"metadata": {},
"source": [
"## Dynamic density correlation\n",
"\n",
"At first order the density correlator contains no interaction vertex, and its deviation from the reference is larger than for $G_0(\\tau)$. The vertex corrections at higher orders bring $\\chi(\\tau)$ progressively closer to the exact result, although it converges more slowly than $G_0(\\tau)$."
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "5fb722e7",
"metadata": {
"execution": {
"iopub.execute_input": "2026-09-29T20:07:12.829925Z",
"iopub.status.busy": "2026-09-29T20:07:12.829785Z",
"iopub.status.idle": "2026-09-29T20:07:12.954077Z",
"shell.execute_reply": "2026-09-29T20:07:12.952336Z"
}
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"for order, r in results.items():\n",
" oplotr(r['chi_tau'], label=label(order))\n",
"oplotr(chi_tau_ed, 'k--', label='pyED')\n",
"\n",
"plt.xlabel(r'$\\tau$'); plt.ylabel(r'$\\chi(\\tau)$'); plt.ylim(bottom=0); plt.legend(loc='best');"
]
},
{
"cell_type": "code",
"execution_count": 10,
"id": "280169e9",
"metadata": {
"execution": {
"iopub.execute_input": "2026-09-29T20:07:12.956946Z",
"iopub.status.busy": "2026-09-29T20:07:12.956740Z",
"iopub.status.idle": "2026-09-29T20:07:12.960940Z",
"shell.execute_reply": "2026-09-29T20:07:12.959620Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"order 1: max|G_0 - G_0^ED| = 3.08E-02, max|chi - chi^ED| = 1.90E-01\n",
"order 2: max|G_0 - G_0^ED| = 1.37E-02, max|chi - chi^ED| = 7.66E-02\n",
"order 3: max|G_0 - G_0^ED| = 1.51E-03, max|chi - chi^ED| = 2.33E-02\n",
"order 4: max|G_0 - G_0^ED| = 3.69E-04, max|chi - chi^ED| = 7.57E-03\n"
]
}
],
"source": [
"for order, r in results.items():\n",
" G_err = np.max(np.abs(r['G_tau'].data - g_tau_ed.data))\n",
" chi_err = np.max(np.abs(r['chi_tau'].data - chi_tau_ed.data))\n",
" print(f'order {order}: max|G_0 - G_0^ED| = {G_err:2.2E}, max|chi - chi^ED| = {chi_err:2.2E}')"
]
},
{
"cell_type": "markdown",
"id": "b5426c8c",
"metadata": {},
"source": [
"The deviation from the exact reference decreases with every expansion order, for $G_0(\\tau)$ faster than for $\\chi(\\tau)$."
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.12.12"
}
},
"nbformat": 4,
"nbformat_minor": 5
}